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Only a few classes of quantum algorithms are known which provide a speed-up over classical 
algorithms. However, these and any new quantum algorithms provide important motivation for 
the development of quantum computers. In this article new quantum algorithms are given which 
are based on quantum state tomography. These include an algorithm for the calculation of several 
quantum mechanical expectation values and an algorithm for the determination of polynomial fac- 
tors. These quantum algorithms are important in their own right. However, it is remarkable that 
these quantum algorithms are immune to a large class of errors. We describe these algorithms and 
provide conditions for immunity. 

PACS numbers: 



I. INTRODUCTION 

There are only a few known quantum computing (QC) 
algorithms which provide a speed-up over their classi- 
cal counterparts. The reasons for this are not completely 
clear However, those algorithms and the associated 

techniques for solving problems efficiently are quite valu- 
able [313) IE]- For example, there are algorithms which 
belong to the same class as Shor's factoring algorithm Q 
which enable the identification of a hidden abelian sub- 
group by using a quantum Fourier transform. There is 
another set of algorithms belonging to the same class as 
Grover's search algorithm which can be applied to a 
wide class of problems where searching a solution set is 
the optimal known problem-solving strategy. Yet another 
class consists of algorithms for simulating quantum sys- 
tems. Simulation algorithms can provide an exponential 
speed-up over any known classical algorithm for a variety 
of quantum systems [1, i, 0, i, i, M, [HI, E E Q ^ 
and are very promising for many applications in the 
physical sciences. These include atomic, molecular, solid 
state, and nuclear simulations and do not necessarily re- 
quire a fully scalable quantum computing device [16|. 

To achieve the speed-ups promised by quantum com- 
puters, a reliable quantum information processing device 
is required. However, noise and imperfections still stand 
in our way. While the active strategies to prevent errors, 
such as quantum error correcting codes [o, [H, [l^, US] 
may, in principle, be universal as claimed, passive pre- 
vention methods have hardware resource advantages. For 
example, decoherencc-frcc subspaces (DFS) and noiseless 
subsystems (NS) [21, 22, 23, 24] are based on the sym- 
metry of the system-bath interaction, so do not require 
active detection and correction of errors. Another passive 



technique, holonomic quantum computation, is robust 
against stochastic errors in the control process [1^ [2^ . 
When conditions are appropriate, passive strategies can 
be applied during the design of quantum algorithms. 

No matter which error prevention strategy is adopted 
however, it is widely believed that entanglement con- 
tributes to the errors in the system and is also the 
resource which is required to achieve the efficiencies 
promised by quantum computers. Here we take advan- 
tage of quantum entanglement in an obvious way in order 
to provide algorithms which solve some problems in poly- 
nomial time on a QC device. Our quantum algorithms, 
some of which are able to calculate quantities which are 
now clearly out of reach for classical computing devices, 
use quantum state tomography (QST) [27| . 

Our QST-based method complements the scattering 
circuit method [H, [29l 30l. Isil. the quantum phase es- 
timation algorithms 32 . |33|. 34 1 it subsumes, and the 
adiabatic method we previously introduced for pairing 
Hamiltonians [l^. As our current method does, these 
other methods scale polynomially in the size of the in- 
put, and this is an exponential speedup over the best 
known classical algorithms for the same task. However, 
our algorithms are exceptional since they are immune to 
a large class of errors and therefore share some features 
with the passive error protection methods of DFS/NS. 
Unlike this previous work, our emphasis is on the error 
resistance of the algorithms. 

Specifically, we first describe an algorithm for the de- 
termination of various types of observables that one may 
want to extract from a system which is being used as a 
quantum simulator. Second, we provide a method for the 
determination of the factors of a large polynomial using 
a quantum computing device. Our objectives are 1) to 
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show why these algorithms are robust against errors and 
2) how these and other algorithms can take advantage of 
such inherent robustness. 

II. ALGORITHM FOR OBTAINING THE 
EXPECTATION VALUES OF AN OBSERVABLE 

Let us first consider the simulation of a quantum sys- 
tem of N subsystems (e.g., particles). It is well known 

that Hamiltonians H of the form H = J2k=i ^k, where 
L is a polynomial in N, and such that efficient quantum 
circuits exist for each term Hk (e.g., when all Hk have 
a tensor product structure or are simply sums of local 
terms), generate unitaries U = exp(ii/t) which can be 
polynomially simulated j35l] . Even random unitary ma- 
trices can be simulated polynomially [s^. Let us denote 
by 0{N^) the simulation cost of such efficiently simu- 
latable unitaries, where fc is a fixed integer. This then 
yields an efficient algorithm for obtaining the quantum 
state IV') = t/|000...0)jY. 

The most general operator of non-identical d-level par- 
ticles (qudits), up to two-body interactions, can be writ- 
ten as 

0= 5] Oa/3,^5(*j)|a,/?j)(7,<5,|, (1) 

where i,j label the N subsystems and a,/?, 7, (5 £ 
{0,l,...,(i— 1} the states of the qudits, and {|ai)} is 
a basis for the Hilbert space of one qudit. For example, 
this could be a Hamiltonian or a unitary gate. 

We are interested in the expectation value (O) = 
{tp\0 IV") in a given quantum state \^) = U |V'(0)) where 
|'(/'(0)) is an initial state. The existing classical algorithms 
for (O) require the simulation of the unitary matrix U 
with X independent elements. Clearly, the classi- 
cal simulation cost grows exponentially with the input. 
An efficient and general quantum method for obtaining 
(O) when O is unitary is the "scattering circuit" [1^: one 
prepares an ancillary qubit in the state (|0) + |l))/\/2, 
interacts the main system with it using a controlled-[/ 
operation, then measures the Pauli operator tr"*" on the 
ancilla; this yields (O) for qubits [11, [H, US [U (we 
are unaware of a generalization of this method to qudits, 
though believe this is possible). The scattering circuit 
method includes quantum phase estimation algorithms 
[3^ m, [11] as special cases; its computational cost is 
0{N^). Here we introduce a different general method, 
based directly on QST [27]. Our method has a computa- 
tional cost that is higher by a factor of 0{N'^) than the 
scattering circuit, but it does not require an ancilla and, 
more importantly, exhibits a remarkable inherent fault 
tolerance to decoherence errors. 

To this end it is convenient to re-express (O) in a 
form relevant to QST. The two-qudit reduced density 
matrix p'^ is given by p^'^ — J2m ('^IV') (V'I'^)j with 
m running over all the c?^~^ orthonormal basis vec- 
tors, excluding qudits i and j. p*-' is x di- 



mensional, with elements p^^g ^/s — ili^jl P^'' I'^if^j) = 
Em il^^j^W {^\mai(3j) = (ijlma.Pj) (7,(5jm]V') = 
(iplaiPj) (jiSjlip), where we have used that {ip\mail3j) 
are c-numbers and J2m 1"^) ("^1 = 1- Using Eq. ([T]) we 
thus have 

(O)- E 0.f,,,s{^J)p%^^Y.^TiO{^J)p^^). (2) 

This expression implies an efficient quantum algorithm 
for (O), as follows: (0) Classically calculate the d^ x d^ 
matrix elements Oap,^5{ij) for all N{N — l)/2 distinct 
pairs of qudits, in the fixed basis {]ai)}. (i) Propagate 
]V'(0)) to ]V') using J7, which can be done in 0{N^) steps 
as noted above, (ii) Using QST find the c?^ — 1 real com- 
ponents of p*^ (for a given pair of qudits (iii) Re- 
peat steps (i) and (ii) for all N(N — l)/2 distinct pairs of 
qudits. (iv) Repeat step (iii) M times, to obtain an es- 
timate of (O) with a precision (standard deviation) that 
scales as l/-\/M(central limit theorem), (iv) Classically 
evaluate Ey Tr(0(ij)p*''). The total simulation cost is 
0{d'^MN^^^). However, we may note that this might be 
improved using more recent QST methods [37]. 

Note that this method can be generalized to the case of 
n-local observables with many-body correlations. Specif- 
ically, any operator on N qudits can be expressed as 
a linear combination of terms, each of which is a ten- 
sor product of N generalized Pauli matrices (e.g., the 
"very nice error operator basis" ^3]), where we in- 
clude the d X d identity as a generalized Pauli matrix. 
If each of these tensor products contains at most n gen- 
eralized Pauli matrices not equal to the identity then 
the operator is said to be n-local. In the case of a n- 
local operator the obvious generalization of Eq. ^ is 
(O) = E^l^,...^„Tr(0(^l^2...^„)p^^^=•■""), where p'l'^-'- 
is the d" X d" dimensional reduced density matrix of 
particles 1112-- -in- Its d^" — 1 real components can be 
obtained via QST, again using a fixed number M of 

copies of ]'0(O)). This must be done for all ^ ^ ^ 

tuples of particles. Therefore the total computational 
cost of our algorithm for the expectation value in the 
case of n-local observables is 0{d'^'^MN^^'^). The mea- 
surement error e = ((O — C'ost)^)ave (where Cest is the 
estimator employed and averaging is with respect to the 
M repetitions) satisfies the generalized uncertainty rela- 
tion (derived from the Cramer- Rao bound) [4Q] : eAH > 
l/(2\/M), where AH = {(H)^ - (ij2))i/2 jg the variance 
of H on the input state ]V'(0)). This bound is indepen- 
dent of n but depends implicitly on d through AH. It is 
important to note that if ]'0(O)) is itself an entangled state 
of P identical copies then the measurement error can be 
reduced by a factor of P (the Heisenberg limit); the de- 
tails and a general proof of optimality of this bound, as 
well as its achievability, are discussed in [ZH ]. 

An important special case is when is an eigen- 
state of a Hamiltonian. The energy spectrum may then 
be found by preparing a (complete) set of eigenstates 



3 



l^/i^) and measuring the set of expectation values (©)„ = 
{ipnlHlipn) = En- Let US commeut on precision issues in 
this context. Our QST-based method complements the 
scattering circuit method [1^, [29', and the quan- 

tum phase estimation algorithms [32, .33,, ,34,] it subsumes, 
and the adiabatic method we previously introduced for 
pairing Hamiltonians fl5'|. As in our current method, 
these other methods scale polynomially in N, and this 
is commonly considered an exponential speedup over the 
currently known best classical algorithms for the same 
task. However, for error e (defined above) the number of 
digits of precision I in the result is r log(l/e), and both 
the scattering circuit and the adiabatic methods require 
poly(l/e) elementary steps to obtain this precision, due 
to the use of the (quantum) Fourier transform at the mea- 
surement ,42i] . In contrast, an efhcient algorithm would 
only require poly(log(l/e)) number of steps. As observed 
by Brown et al. [42| . while this has no impact for fixed 
precision, the 1/e scaling does imply an exponential scal- 
ing with the number of digits of precision. The origin of 
the /~ log(l/e) scaling is illucidated by Giovannetti et al. 
[4l| , who show that this scaling cannot be improved even 
using entanglement. Namely, they show that entangled 
measurements do not help, and the use of P entangled 
input probes gives at most the Heisenberg limit e~l/P, 
and on the other hand F logP. Thus r log(l/e). While 
our QST based method does not employ a Fourier trans- 
form at the measurement, the general arguments used in 
(4ll | apply to QST as well, so that our present algorithm 
does not improve on the precision issue. As discussed in 
p3 |. the origin of the poly(l/e) number of steps is in the 
use of the Trotter formula for the simulation of U. Use 
of the Solovay-Kitaev theorem [which improves the Trot- 
ter poly (1/e) scaling to 0(log (1/e)) scahng] does not 
help when a fault tolerant implementation is considered, 
since the latter once again leads to the poly (1/e) scaling 
However, it is important to note that these general 
bounds do not preclude specific Hamiltonians from being 
efficiently simulatable in terms of precision requirements; 
indeed the exponential precision slow-down is avoided in 
Shor's algorithm due to the manner in which modular 
exponentiation is carried out Q. Another observation 
is that in some cases it is possible to prepare the final 
state lip) by means other than quantum simulation, e.g., 
via cooling to the ground state, or via adiabatic evolu- 
tion. There are certainly examples where then reaching 
1^) from \^{0)) requires poly(log(l/e)) steps. 



III. ALGORITHM FOR OBTAINING THE 
EXPECTATION VALUES OF AN OBSERVABLE 
OF A FERMIONIC SYSTEM 

It follows from the Jordan- Wigner transformation [i^l 
that there is one-to-one correspondence between fermions 
characterized by the fermionic creation and annihila- 
tion operations and Cj, where j denotes a fermionic 



mode, and qubits: c] (~1)^ ^ (®/=i "^f ) '^t 

[a^ — [a^ ± iay)/2 and cr^'i''^ are the Pauli matri- 
ces]. Therefore, the above algorithm can be applied 
to a fermionic system. However, a one-body or two- 
body interactions of fermions usually corresponds to a 
many-body interaction of qubits, for instance c\cj <^ 

(^izl fjf ^ af ...a^^_^a^ aj where j > i. Even so, 

it is clear that obtaining the expectation value of an ob- 
servable with many-body correlation still requires only 
polynomial time. For example, measuring an observable 
such as (jf(J^(jfj^i...(j^j_i can be accomplished efficiently, 
as long as the distance between i and j is finite and inde- 
pendent of A^. In some cases partial QST suffices to ob- 
tain a desired expectation value. For example, consider 
the one-body Fermi operator h = '^eiUi [iii — c|ci), 
where we assume the e^ are known. Then 

/(V'l/ilV')/ = ^e^/ (V'lni IV'); 

= E^^/'li' (3) 

where = E^Ii Ea^^o I") (" = {ai, aAr}) 

is an arbitrary pure state of N qubits and — 
(Tr^-^j ]■!/)) (■01) 11, j = l...iV. <^ means that there is a one- 
to-one correspondence between f expressed by Fermi 
creation operators on the vacuum state and expressed 
by the superposition of computational bases [44| . 



IV. INHERENT ROBUSTNESS 

The accuracy of a usual quantum algorithm requires 
that the final wave function or density matrix be 
Pi — iV'o) ("Vol- In reality, due to errors, the actual density 

matrix will be given by pA = ^k=oPk I'/'fe) ("^fcl which is 
different from the ideal one pj. However as long as the 
following relations are satisfied; 

P' = Trj/5/ = Tr^p^ or, p'^ TiijPi Tr^^p^ (4) 

our algorithms will give the same results, where the sub- 
scripts i and j means trace over all degrees of freedom 
excluding i and j. There are only 3A'^ or 15N'^ constraints 
respectively. This implies that our algorithms are much 
more fault-tolerant than a generic one. The reason that 
the algorithms are more robust is that there are 2^ in- 
dependent coefficients in pA- However, we only require 
that the above relation holds independent of the other 
various parameters in the system. 

As a motivational example, suppose an expected final 
state is IV'o) = a [00) + 6 Jll) , but due to dephasing errors, 
we actually get 

PA = lal2 joo) (OOl + JfoJ^ 111) (111+C 100) (lll+C* 111) (GO] 
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where C is an arbitrary number and is zero when com- 
plete phase damping occurs. No matter what value of 



|6p 



Therefore, dephasing does not affect the validity of our 
algorithm in this case. The algorithms are, to some ex- 
tent, self-protected. 

To provide some general conditions under which our 
algorithms are robust, let us start with some definitions. 
Let 



PA<i$ PB <i$ PE, 



(5) 



where pA is the subsystem we wish to study, ps is the 
rest of our system, and ps is the density operator for the 
environment. We can assume that each of these is a pure 
state and the whole system plus environment is pure and 
initially completely separable. Now. let [/ (Xi be the 
ideal unitary operation for our simulation algorithm and 



p/ = [/ IepU^ O Ie 



(6) 



Let V be the non-ideal operation. We can write the 
condition for the algorithm to give the same result for 
the expectation value of an operator O as Tr{Op') = 
Ty{OWp'W'<), where W = V{W ® Ie)- Let the basis for 
the algebra of operators be traceless and Hermitian and 



represented by X. 



1,2,3, for subsystems A,B,E 



respectively with a E {1, ...,<P — 1}. Then we may write 
the density operator for the ABE system as 



PABE = IaBE + + ^"^^IZ^^^i ^ 



+Iab (8) y^CfcA 



(3) 
k ■ 



(7) 



If O = n • A(i) (or C = n • A^^) O Ie), then O acts as a 
projector onto the subspace A and the expectation value 
of O is 



(O) = TiiOWp'W) = Tr(n • X^^^Wp'W). 



(8) 



So if Wp'W^ has the form pabe, then {O) = n-a, where 
a = {ai,a2, . . . }. Likewise, if 



/ = Iab + Y1 «i \ ®Ib + Ia(^Y1 ^'j^. 



/\(2) 



(9) 



completely positive (CP) maps. Let us choose an initial 
density matrix p which will be acted upon by a CP map 
corresponding to the operator-sum decomposition with 
operators Ai. We then want to find: 



(0)i=Tr(0^A,p4). 

i 

Note that this can be written as 

{0),=Tr{J2AlOAip) 



(10) 



(11) 



so that the condition for the same result to be obtained 
from a different set of operators Bi is 

= - {0)2 

= Tr ^AjOAip^ -Tr (^B}OB,f}j 

Y,A\OA,-Y,B]OB,\p 



= Tr 



(12) 



Therefore, we may also say that the expectation value is 
invariant under transformations which arc comprised of 
the little group of O. This is true for both the unitary 
description above, as well as the operator-sum decompo- 
sition. 

Let us simplify to the case of a qubit. Letting O be 
traceless and Hermitian and Bj = f3jl + bj • a and Ai = 
ail + Si - a we may obtain the relation 

Y^AlOA, = Y^^ilidi xa*)-n 

i i 

-|-(|ai|^ - a ■ a*)nt(Tt 

+ [ai{a* X n)t - a*{di x n)t](7t 

+ [{d* ■ n)ait + {Si ■ n)a*^]at), (13) 

where the sum over t is implied. Simplifying further by 
letting O = CT3 and p = (l/2)(/ -|- a^), we can write the 
condition as 

^[i(afc X 0^)3 + {\akf - \akif - |afc2|^ + lafcsl^)] 

fe 

-Y^m X b*)s - - \bn\' - \bj2\' + \bjsn] = o. 



then Tr(Op') = n-a^ . Therefore, for these to be equal, we 
require that n-a = n-a^ .liwe write n-a^ = |n||a^| cos9^ , 
andn-a= |n||a| cos ^, then we need |a''^| cos 0^ = \a\ cosO. 
For a two-state subsystem A, this leaves one degree of 
freedom, the little group of the vector a. This is stated 
in terms of the coherence vector for a general expectation 
value for a rf-state system. 

We may also show that the robustness can be expressed 
in terms of the expectation value of the operator O and 



Note that J^i^l^i = ^. implies Eid^iP + ' «i) = 1> 
and J^il^i^it +ctiait + i{ai x a|)t] = 0. So the result can 
be expressed in terms of two equations 

Y^[i{ak X al)3 - 2(|afeip + \ak2\^)] 

k 

-Ymxb*)s + 2{\bi^f + \bj2f)]=0. 



5 



and 

k 

j 

To summarize, owr simulation algorithms, based on quan- 
tum state tomography and the expectation value of an op- 
erator, are immune to errors which act as the little group 
of transformations of the initial density operator or the 
operator for which we seek the expectation value. 

V. ALGORITHM FOR FACTORING A 
POLYNOMIAL 

We now present one more algorithm which can be 
implemented via state tomography and which is robust 
against the aforementioned class of errors. Consider vari- 
ables Xi,yi , where i = 1, 2, iV, and a class of homoge- 
neous functions spanned by the set of products of Xi,yi. 
For instance when N — 2, the set is xiX2, Xiy2, yiX2 and 
2/12/2- There is a one-to-one correspondence between this 
set and the computational basis for two qubits. The 
linear combination of the set defines a class of homo- 
geneous functions. For instance, consider the two func- 
tions X1X2 + yi2/2 and xij/2 + J/i?/2- The former cannot be 
factored, while the latter can be factored into the form 
(xi -|- yi ) X ?/2 • In some circumstances, it may be easy to 
tell whether or not this can be factored, if we know the 
concrete form of the homogeneous function. However, if 
a homogeneous function contains many terms, in general 
it will become difficult. Consider such a function derived 
from a matrix U acting on a basis set such as X1X2X3...XN 

fN{xi,yi) ^ UX1X2X3...XN (14) 

where 17 is a 2^ x 2^ matrix. To represent the function, 
a classical computer needs to handle 2 x 2^ x 2^ inde- 
pendent numbers in 2-^^ x 2^ complex matrix elements of 
U in order to simulate it. We may assume U is unitary 
so that it preserves the norm of the function. However, 
this still requires a classical simulation of 2^ x 2^ in- 
dependent numbers in the matrix U . When N = 300, 
approximately 10^^° independent numbers must be han- 
dled. 

Given a unitary matrix U which could be a random 
unitary matrix, or the quantum state |/) — U |000...0)^, 
which can be simulated polynomially, say N'' , [36j where 

is a fixed number, we will determine the presence of 
a factor axi -\- hyi. We first obtain the reduced density 
matrix of the ith qubit = (1/2) (/ + n' • c?*) by using 



quantum state tomography, which requires need a fixed 
number, M, copies of |/), as discussed above. Then, we 
calculate the von Neumann entropy of p*. If the entropy is 
zero, the ith qubit is separable from the others, meaning 
that there is a factor aXi -\- hyi in fN{xi,yi) where a — 
^Jl + nl, and b — {n\.-\-in\^) / 1 -\- given by the matrix 
elements of p*. Otherwise, there is no such factor. The 
total number of steps in the quantum procedure is MN^ . 
The same procedure can be used to find higher order 
factors. For example, a, factor aXiXj-\-bxiyj-\-cyiXj-\-dyiyj 
can be found by measuring the reduced density matrix 
p^^ for the ith and jth qubits. 

Furthermore, the generalization to many-qudit sys- 
tems, (each subsystem has an arbitrary dimension), and 
thus multivariate polynomials can be accomplished by 
using the generalized coherence vector, or generalized 
Bloch vector [H, IH, li^l . Let A[ (g) ^| be a Hermitian 
basis for a system of coupled qudits, with arbitrary di- 
mensions for all components. Let form a basis for the 
ith subsystem with Tr(A[A*) = 2Srt and fj,j a basis for 
the rest of the system. Then, given the state C/|000...0) 
for the whole system, the corresponding reduced density 
matrix for the ith subsystem has p — (l/d)(l + m • A). 
The conditions for the system to be factorisable with re- 
spect to the ith subsystem is that m • m = N{N — l)/2 
and dY^mrmsN{N - 1)/{2N - 4) = wt, with = 
(7V/4)Tr({A[,Af}A*). As usual, {•,•} denotes the anti- 
commutator. These conditions indicate that the qudit is 
in a pure state and thus has zero entropy. Using the Her- 
mitian basis for the operators in this protocol provides 
an explicit measurement basis for the identification of the 
reduced density matrices. 



VI. CONCLUSIONS 

In this paper, we introduced quantum algorithms 
based on quantum state tomography. The simulation 
algorithms are clearly polynomial while the best known 
classical counterparts of the simulation algorithms are 
exponential. We suspect that the polynomial factoring 
algorithm is also more efficient although we have not 
proved this generally. Certainly in the case that the 
unitary U must be simulated, we achieve an exponen- 
tial speed-up. We emphasize that the algorithms are, to 
a large degree, self-protected against a large class of er- 
rors. This work brings together two important aspects 
in quantum information science, algorithms and passive 
correction. We expect that the family of quantum algo- 
rithms which are error-avoiding algorithms should receive 
much more attention in future studies of quantum algo- 
rithms. 
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